clear all
close all
[data params]=fca_readfcs('C:\Documents and Settings\Michael Ferry\My Documents\experiments\MF216\MF216\MFSC101_12.fcs');
indicesToProcess=11;
randIndices=randint(1000,1,[1 10^5]);
smooth_value=11;
smoother='sgolay';
%hist_res=2^8;
%numbins=params.par(indicesToProcess).range;
%scaled_bins=exp((1:hist_res)*(log(numbins)/hist_res));



%hist(data,scaled_bins); 
% [x nX]=hist(data(:,indicesToProcess),2^12);
% % figure
% % bar(nX,x);
% figure
% [indX cX]=kmeans(data(:,indicesToProcess),2,'dist','sqeuclidean');
% %[silh2,h] = silhouette(data(randIndices,indicesToProcess),indX(randIndices),'sqeuclidean');
% x1Lower=min(data(find(indX==1),indicesToProcess));
% x1Upper=max(data(find(indX==1),indicesToProcess));
% x2Lower=min(data(find(indX==2),indicesToProcess));
% x2Upper=max(data(find(indX==2),indicesToProcess));
% figure
% bar(nX,x);
% line([x1Lower x1Lower],[0 1000],'Color','r');
% line([x1Upper x1Upper],[0 1000],'Color','r');
% line([x2Lower x2Lower],[0 1000],'Color','g');
% line([x2Upper x2Upper],[0 1000],'Color','g');

% figure
% hist(data(find(indX==1),indicesToProcess),2^8);
% figure
% hist(data(find(indX==2),indicesToProcess),2^8);

%[y1 nY]=hist(log(data(:,indicesToProcess)),2^8);
[y nY]=hist((data(:,indicesToProcess)),2^12);
scaler=log(1:2^8);
%y=y1./scaler;
r=smooth(y,smooth_value,smoother);
figure
%bar(nY,r)
hist(log(data(:,indicesToProcess)),2^8);
sensitivity=85;
peaks =fpeak(nY,r,sensitivity,[2 max(nY)-1 10 Inf]);
for i=1:length(peaks(:,1))
     %line([peaks(i,1) peaks(i,1)],[0 200],'Color','r');
end
%numOfpeaks=length(peaks(:,1));
numOfpeaks=2;
[indY cY]=kmeans(log(data(:,indicesToProcess)),numOfpeaks,'dist','sqeuclidean');
%[silh2,h] = silhouette(log(data(randIndices,indicesToProcess)),indY(randIndices),'sqeuclidean');

%y1Lower=log(min(data(find(indY==1),indicesToProcess)));
%y1Upper=log(max(data(find(indY==1),indicesToProcess)));
%y2Lower=log(min(data(find(indY==2),indicesToProcess)));
%y2Upper=log(max(data(find(indY==2),indicesToProcess)));
% figure
% bar(nY,y);
for i=1:numOfpeaks
    yLower=log(min(data(find(indY==i),indicesToProcess)));
    yUpper=log(max(data(find(indY==i),indicesToProcess)));
    line([yLower yLower],[0 1000],'Color','g');
    line([yUpper yUpper],[0 1000],'Color','g');
end


% figure
% bar(nX,x)
% line([exp(y1Lower) exp(y1Lower)],[0 1000],'Color','r');
% line([exp(y1Upper) exp(y1Upper)],[0 1000],'Color','r');
% line([exp(y2Lower) exp(y2Lower)],[0 1000],'Color','g');
% line([exp(y2Upper) exp(y2Upper)],[0 1000],'Color','g');

%[indZ cZ]=kmeans(data(:,2),2,'dist','sqeuclidean');
% title('histogram of raw data')
% y=zeros(length(x),1);
% y(1:2^6)=smooth(x(1:2^6),31,'sgolay');
% y(2^6+1:2^8)=smooth(x(2^6+1:2^8),51,'sgolay');
% y(2^8+1:2^12)=smooth(x(2^8+1:2^12),101,'sgolay');
% %y(10001:2^12)=smooth(x(10001:2^12),151,'sgolay');
% figure
% bar(y)
% title('smoothed histogram of x')
% 
% z=diff(y);
% figure
% bar(z);
% title('differential of smoothed x')
% %z=smooth(y,11);
% % figure
% % bar(z);
% 
% sign_prev=sign(z(1));
% cluster=1;
% ind=zeros(length(z),1);
% for i=1:length(z)
%     if(sign_prev==sign(z(i)))
%        ind(i)=cluster; 
%     else
%         if(abs(z(i))>10)
%             sign_prev=sign(z(i));
%             cluster=cluster+1;
%             ind(i)=cluster;
%         else
%             ind(i)=0;
%         end
%     end
%     
% end
% cluster        
    